Домашняя работа по ТВИМС

Учебный ноутбук с R(t) моделями

Выполнил:

Вольхин Данил Федорович

Email:

dfvolkhin@edu.hse.ru

Дата:

30 октября 2024

😊 Приятного просмотра 😊

Учебник по генеративным моделям

В этой лекции вы узнаете:

В этом блокноте объясняется модель, на которой работает https://rt.live с ее прогнозом эффективного числа воспроизводства. Она основана на версии 1.0.0 модели.

Что такое Rt?

$R_0$ («R-ноль») описывает фактор воспроизводства заболевания — т. е. скольким другим людям один инфицированный человек передает болезнь. Если это количество больше 1, то у нас на руках эпидемия, как в случае с COVID-19.

Однако $R_0$ предполагает, что не предпринимаются никакие контрмеры для сдерживания распространения вируса. Таким образом, более важной мерой для отслеживания является $R_e(t)$ — изменяющийся во времени эффективный фактор воспроизводства, т. е. сколько людей заражает один человек в определенный день $t$.

По мере введения мер изоляции и социального дистанцирования мы ожидаем, что это количество снизится, в идеале ниже критического количества 1, потому что тогда со временем болезнь просто сойдет на нет.

Обычно мы извлекаем $R_e(t)$ из чего-то вроде модели SIR (восприимчивый-инфицированный-выздоровевший) или SEIR (которая добавляет категорию подверженных воздействию), которые являются классическими эпидемиологическими моделями отсеков.

Модель SIR — это то, что rt.live использовал в начале. Однако модели SIR/SEIR также являются лишь приближениями к реальной вещи и содержат довольно много встроенных предположений. Текущая модель проще и делает меньше предположений. Кроме того, модель SIR описывается как ODE, что вызывает различные технические проблемы. Решение ODE требует довольно много времени, и хотя существуют приближения в замкнутой форме, которые работают быстрее, мы обнаружили, что они довольно ненадежны.

Вместо этого текущая модель использует простую генеративную логику для объяснения того, как начальный пул инфицированных людей распространяет болезнь в каждый момент времени в соответствии с текущим фактором воспроизводства.

Генеративная модель

Нам нужен пакет из репозитория rtlive для загрузки данных и получения наших дистрибутивов.

Импорт библиотек

Задаем параметры

Загрузка данных

Хотя раньше мы рассматривали кумулятивные подтвержденные случаи, имеет смысл рассматривать новые ежедневные случаи. Это также обходит проблему, что все наши предыдущие модели допускали снижение случаев со временем, в то время как процесс генерации данных не допускал этого.

Предположим, что для идеальной болезни мы начинаем с одного инфицированного пациента (первичная инфекция) в день $0$, который заразен только один день и в этот день заражает 2 человек (вторичная инфекция), которые заболевают на следующий день. Таким образом, эта болезнь имеет фактор воспроизводства $R_0$, равный 2. Мы могли бы записать, что в день $t$ число вновь инфицированных $y_t$ равно:

$$ y_t = y_{t-1} \cdot R_0 $$

Довольно просто. Эта логика порождает классическую модель экспоненциального роста, которую мы видим в эпидемиях:

Построение модели

Чтобы встроить эту рекурсивную природу в PyMC3, нам нужно создать цикл внутри нашей модели. Поскольку PyMC3 компилирует все в theano, мы не можем просто использовать цикл Python. К счастью, theano имеет собственную логику цикла в форме scan. К сожалению, синтаксис довольно запутанный. Не беспокойтесь слишком сильно обо всех отдельных аргументах, но обратите внимание, что наша генеративная логика реализована как лямбда-функция, переданная в fn. scan затем будет проходить по sequences и обновлять y на каждой итерации.

Здесь мы также впервые используем Deterministic. Все, что он делает, это сохраняет переменную вашей модели в трассировке. Если у вас есть какое-либо выражение, состоящее из других случайных величин, часто удобно записывать их вместе с апостериорными параметрами.

Изменяющееся во времени воспроизведение

Однако, как мы уже обсуждали, нас больше интересует эффективная скорость воспроизведения как функция времени $R_e(t)$. Мы можем просто переключить это в нашу идеализированную генеративную модель:

$$ y_t = y_{t-1} \cdot R_e(t) $$

Вы можете видеть, что в последний день, когда $R_e(t)$ равен 1, мы получаем то же количество новых заражений, что и в предыдущий день, потому что каждый зараженный заражает всего одного человека.

Неявное предположение в генеративном процессе выше заключается в том, что зараженный человек заразен только один день, и что затем ему требуется всего один день, чтобы заразить других людей.

Построение модели

Какой тип априорной вероятности будет установлен для этого изменяющегося во времени фактора воспроизводства $R_e(t)$? Априорная информация, которая у нас есть, заключается в том, что она не будет резко меняться от одного дня к другому. Мы кодируем эту информацию с помощью априорной вероятности случайного блуждания:

$$ R_e(t) \sim \mathcal{N}(R_e(t-1), \sigma) $$

где $\sigma$ контролирует ширину шага случайного блуждания. К счастью, это распределение уже включено в PyMC3. Поскольку мы знаем, что фактор воспроизводства может быть только положительным, мы берем exp неограниченного случайного блуждания (которое тогда будет в логарифмическом пространстве), чтобы сделать его положительным.

Таким образом, модель становится:

Учет задержки заражения

В действительности время, необходимое первичному человеку для заражения других, следует распределению. Они могут заразить одного человека на следующий день, двух — через день и т. д. Это распределение задержки официально известно как «время генерации», и мы смоделируем его с помощью распределения вероятностей из этого исследования30119-3/pdf) (исследование на самом деле дает оценку чего-то, тесно связанного со временем генерации, но ее можно использовать для этого).

В статье приведена простая формулировка того, как вывести распределение времени генерации для COVID-19:

Чтобы включить этот эффект в нашу генеративную модель, нам нужно сделать свертку. Интуитивно понятно, что вместо новых случаев в день $t$, зависящих только от новых случаев в день $t-1$, теперь они зависят от новых случаев в (потенциально) все предыдущие дни, поскольку между моментом заражения человека и заражением другого человека могло пройти 5 дней. Нам нужно учесть всех этих ранее инфицированных людей и то, с какой вероятностью они заражают людей сегодня. Мы достигаем этого, взвешивая количество новых инфицированных людей $i$ дней назад -- $y_{t-i}$ -- по времени генерации $g_i$ для этой конкретной задержки, а также по эффективному числу воспроизводства в этот день $R_e(t-i)$:

$$ y_t = \sum_{i=1}^{M}y_{t-i} R_e(t-i) g_i $$

Более подробную информацию об этом генеративном процессе см. в этой публикации: https://staff.math.su.se/hoehle/blog/2020/04/15/effectiveR0.html.

Обновив нашу модель генеративного процесса соответствующим образом, мы получаем:

Как видите, учет задержки между передачей болезни от одного человека другому значительно замедляет распространение. Чем больше время генерации, тем медленнее распространение.

К сожалению, здесь все становится немного техничным, не беспокойтесь слишком сильно об индивидуальном коде, но поверьте, что мы включаем время генерации сейчас в нашу рекурсивную функцию scan.

Апостериорный прогноз

Получение числа инфицированных

Пока что у нас есть генеративная модель того, как люди передают болезнь от одного человека к другому. Однако у нас нет данных о том, когда люди передали болезнь, у нас есть данные о том, кто получил положительный результат теста. Поэтому нам нужно еще больше отсрочить эту функцию, когда инфицированный пациент фактически появится в качестве положительного теста в наших данных.

Для этого мы будем использовать распределение задержки между заражением и подтвержденным положительным тестом, также известное как распределение задержки начала. Чтобы оценить это распределение, мы можем использовать данные из Рабочей группы по открытым данным COVID-19, которая спрашивала пациентов с COVID-19, как давно у них начались симптомы (Xu et al., 2020). Однако начало симптомов не то же самое, что и при заражении, потому что есть инкубационный период, в течение которого вирус распространяется в организме, пока симптомы не заметны. Чтобы исправить это, мы просто добавляем инкубационный период к началу распределения задержки начала (rt.live предполагает 5 дней). Вы можете увидеть это смещение в плоской области с 1 по 5 день на следующем графике.

Теперь нам нужно лишь свернуть полученную выше функцию количества людей, заразившихся каждый день, с распределением задержки начала заболевания, чтобы получить количество людей, которые обнаружат положительный результат теста в этот день.

К сожалению, это не совсем сработало. Я не уверен, что не так с моделью, но я оставил это здесь, чтобы продемонстрировать, какие вещи можно делать внутри модели. Исправление модели оставлено в качестве упражнения для мотивированного читателя ;).

Корректировка количества проведенных тестов

При рассмотрении количества сырых положительных тестов становится ясно, что на это число будет влиять количество людей, которых вы протестировали: чем больше вы тестируете, тем больше случаев вы обнаружите.

Это важно для моделирования, поскольку существует огромная изменчивость в количестве проводимых тестов с течением времени (наращивание тестов в целом по мере создания большего количества возможностей для тестирования, но также и потому, что по выходным обычно проводится меньше тестов). Это сместит нашу оценку $R_e(t)$, если не учитывать это.

Таким образом, в модели мы умножаем воздействие теста $e_t$ (нормализованную величину, пропорциональную количеству проведенных тестов) на количество положительных тестов из генеративного процесса. Интуитивно, если мы тестируем вдвое больше, мы ожидаем, что положительных тестов будет вдвое больше.

Таким образом, ожидаемое количество положительных тестов $\tilde{z_t}$ будет:

$$ \tilde{z_t} = z_t \cdot e_t $$

где $z_t$ — это выход генеративной модели с примененными задержками.

Подводя итог генеративному процессу

  1. Происходит первичное заражение (это момент времени, к которому мы хотим отнести Rt).

  2. Время генерации проходит до тех пор, пока не произойдет вторичное заражение.

  3. Время начала проходит до тех пор, пока у вторично инфицированного человека не появятся симптомы и тесты не станут положительными. Это количество положительных тестов, которое мы ожидали бы, если бы тестирование было постоянным.
  4. Умножьте количество положительных тестов (если бы тесты были постоянными) на воздействие тестирования, чтобы получить количество ожидаемых положительных результатов. Это выход модели, который мы используем для подгонки данных.

Корректировка прогнозов инфицированных с учетом вероятности положительного теста

test_adjusted_positive = pm.Deterministic(
    "test_adjusted_positive", conv(infections, _p_delay, len_observed)
)
# Определяем случайную величину по распределению Дирихле с alpha=[1,1,1,1,1,1,1]
w = pm.Dirichlet("w", a=[1, 1, 1, 1, 1, 1, 1])

# Применяем `scan` для итерации по `test_adjusted_positive`
weighted_output, _ = scan(
    fn=lambda t, test_value: test_value * w[t % 7],                   # Функция для итерации
    sequences=[tt.arange(test_adjusted_positive.shape[0]), test_adjusted_positive],  # Временные шаги и данные
    non_sequences=w                            # Параметры, не зависящие от шага (веса w)
)

# Сохраняем результат в детерминированной переменной для анализа
weighted_infections = pm.Deterministic("weighted_infections", weighted_output)

Здесь нужно просто домножить test_adjusted_positive на w [t%7], где w это Дирихле

Интерактивное использование модели

Давайте посмотрим, как на самом деле использовать модель интерактивно. В качестве примера региона мы используем Массачусетс.

Чтобы запустить модель, мы просто создаем ее экземпляр и вызываем .sample():

summarize_inference_data() переформатирует необработанную трассировку выборки, чтобы с ней было легче работать.

График выше показывает, что именно происходит на каждом этапе. Давайте применим логику в обратном порядке, на этот раз перейдя от данных к скрытым причинам. Во-первых, у нас есть наблюдаемые данные «Сообщенных положительных тестов» серого цвета. Это то, что наша модель пытается объяснить с помощью «Ожидаемых положительных тестов», выходных данных модели, которые выводятся из «Ожидаемого числа положительных тестов, если бы тесты были постоянными». Это число мы сворачиваем с распределением задержки начала, чтобы получить из «Начала инфекции», т. е. когда люди фактически заразились. Оттуда мы используем распределение времени генерации, чтобы прийти к тому, когда произошла фактическая передача заболевания.

Наша оценка $R_e(t)$ затем основана на том, когда произошла фактическая передача, а не на том, когда мы наблюдали положительный тест: